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SOME PROBLEMS OF HEAT TRANSFER IN LIQUID-COOLED REACTORS 

V.S.Os machhin 
1 . INTRODUCTION 

r lhe modern nuclear power reactors are highly forced. Large heat fluxes and specific power, 
rigid limitations on the coolant and fuel elements temperatures demand knowledge of tempera- 
ture distribution in the core with a high reliability. 

The main problem in thermal reactor design is to prove reliably that the impermissible 
coolant and fuel element temperatures will not be realizedat any reactor operational regimes. 

Specific features of heat transfer problems in modem reactor technology are defined with 
core arrangement, which for the most of power liquid cooled reactors is a set of assemblies 
of canned cylindrical fuel elements. The knowledge of the flow pattern and heat transfer rate 
is necessary for determination of the temperature distribution in such a complicated geometry. 

In compact cores it needs to account for the effect of the axial heat flux variations on the 
heat transfer coefficient. These problems are considered in the paper. 

2 . TRANSPORT PROCESSES IN TURBULENT FLOW 

The motion of fluids in nuclear reactors is turbulent. Theoretical study of the processes 
occurring in turbulent flows is very complicated/ Due to irregularity of turbulent motions it is 
desirable to applicate statistical methods. However, the realization of such an approach in 
all details is very difficult. Therefore, in practice semi-empirical phenomenological theories 
are used. The object of these theories is to find some relations between mean and fluctuating 
parts of the motion. On PrandtPs mixing-length theory the velocity pulsations in a flat channel 
can be presented as 


Ul ~ ly c!y ’ 


where 1. 


1 y - is a characteristic distance passed by the pulsations. 

Let us derive a more exact formulation of PrandtPs hypothesis. Introducing the suitable 
Green functions, we can write the equations of motion for the turbulent fluctuations 


duf 

1 


u} dix} <HJi 1 .1 #ul , 

7 +Uk ^r +u ^^r = _ "^: +l ' — 2 + 1 — ^4 - “j 1 u 1 ) 

C/Xj dxfc p dXj ^9 dxjj * k 1 k 


( 1 ) 
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I hat Green’s function is determined by the following equation 
dV o , 

— +(UV)V 0 -„AV 0 = -SCF-r ] )S(t _ |1) (2) 

with the corresponding boundary conditions. 

The function V 0 describes the velocity distribution in liquid flow under an instantaneous 
local disturbance. 

Then equations (1) can be written as 

u}(r,t) = /v 0 (T, + J — 1 _i- ( ujuJ - uJ u[)] dF^dt 1 (3) 

1 p dX[ * dxfc <?x^ 1 K 1 k 

Such a form of Eq (1) permits an iteration method to be applied. To obtain the first approxima- 
tion formulas we may neglect in Eq. (3) all nonlinear terms and the pressure pulsations. 

Then 

u j (7> 0 — .W 1 — - V Q (7.7 1 , t, t 1 ) dr 1 dt 1 (4) 

Provided the width of the region, where the probability V Q is not zero, is small compared with 
that one, where significant change of the mean velocity occurs, Eq. (4) can be simplified. 
Expanding the mean velocity gradient into a Taylor series near the point 7 and limiting to two 
series terras, we obtain the approximate formulae for the velocity fluctuations 

U JCF, “ f“l (T 1 ( t 1 )V 0 (7,rit. t 1 )dT 1 dt 1 +^- /I7 -7 ^(7 1 ,t 1 )V 0 (7, r^t.djdF 1 

(5) 

In the turbulent core near the flow symmetry axis the first mean velocity derivatives are small. 
Therefore the velocity pulsational components are determined by the second items in Eq.(5), i.e. 
are dependent upon the second mean velocity derivatives. Conversely, far from the flow sym- 
metry axis the first items of the equation are large. In this region 

u i (7 ,0 - /Uj 1 (7 1 , t 1 ) V 0 (7,7 1 I t t t 1 )drllt 1 -l k (6) 


Note, that the quantity 


l( k ' ° = / u k < 7 V o ( r,T 1 , t, t 1 ) dF 1 dt 1 


has the length dimension and represents the width of the region which is a velocity pulsation 
"supplier”. 

To improve Green’s function we must take into account the turbulent diffusion of the 
velocity fluctuations. Using these relations and determining the pulsational diffusion coefficient 
by means of the equality 
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I) = 


u i u k " u i u k 


( 7 ) 


<?*k 


we may introduce the new Green function by means of the equation 


|X + (UV ) v -VDW - V. V = -S ( F -7 1 ) S (t - t 1 ) 
<9t ^ 


( 8 ) 


Then the second approximation velocity pulsation formulae are 

u j’ (f" » 0 -' / • uJ (7 1 ,t 1 )V(r,7 1 ,t,t 1 )dr I dt 1 

i dx k k 

As above, the velosity fluctuations in the region not far from the wall equal 


u k ( r .0=-j^- /uj^r 1 ,t 1 )V(r, r 1 ,t,t 1 )dr W 


d U; 


(9) 


( 10 ) 


Using this expression the turbulent stress tensor components can be calculated 


“M =lira 27^ /u|(7,t)dt/u 1 k (7 1 ,t 1 )V(7 1 7 1 1 t, tbdrGt 1 


( 11 ) 


V-*OQ 


Since Green’s function is dependent on the difference r = t — t ^ only, relation (11) can be 
written in the more convenient form 


.1 1 r 
dx k 

V -► 00 


u! = lim — - /dr — / uj (r . t) dt /u[(r G-r) V(r,r ^rldr 1 = 
1 J d x L . 6 2v-v J k 


d U 

<?x k o 


— / dr/ u!(r ,t)u! (r \t - r) V(r,r\ r)dr^ 


dUj oe *t d U: 

-/K jk (r,„d„J t ._i 


( 12 ) 


T ,- N 


The factor-v.^fr ) •-» / Kj k (r , r)dr is the component of the eddy diffusivity tensor. In the 
Lagrange coordinate system the function v(rj,r) approximately equals v *e ^ r , where 


k = 


a (v + D) 


, a is constant, e is a characteristic dimension of the region. 
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v 


Then for homogeneous isotropic turbulence may be written as the eddy diffusivity of 


momentum 


= f u l(t)ul(t — r)e“k r dr 


This expression correlates with Taylor's determination of the turbulent diffusion coefficient 

[ 1 , 2 ]. 

Similarly using the equations for the temperature fluctuations 

dT 1 ,, ,rrl i <?T ^T 1 s 

- 77 -+^- — +u t^— =a T + H T k 1 ' 

"■ K dXL Jx4 5X1, K 


>x k 8x4 5x k 


and introducing temperature Green’s function by means of the equation 
dW G i i 

? + (UV)W 0 — a A W Q = — $ (7 — 7 *) S ( t - t 1 ) (1! 

d t 

we can obtain the integral equations for the temperature fluctuations in the following form 

T 4 (F ,t) = /W (7 , r“ 1 ,t, t 1 ) [u|~^— —(u^-u^ThdF 1 ^ 1 (ll 

k 5x k * x k k k 

In the first approximation 

T 1 (F ,t) - j'u£ il V^(F,r 1 1 t 1 t 1 )dF 1 dt 1 (1 


or with the same assumptions, which have been used in deriving cf formula ( 5 ) 

thrfl.IL /»i(ri.,i )t v (r . rVjdr 1 * 1 + 




In the region not far from the wall 




T — • 

where 1 ^ ( F, t) is an analog of the mixing-length for the temperature pulsations. Just as in 

equality (7) it can be approximately accounted for the turbulent diffusion of the temperature 
pulsations by means of introducing the turbulent diffusion coefficient 
326 
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Then new temperature Green's function is detennined by the equation 
v + (U V ) W - vn T V W - a AW = - S ( F - F *) S (t - t 1 ) 

at 

Then the temperature fluctuations 

T 1 (F , t) =r /ul ( F 1 , t 1 ) W ( F , F ^t, t 1 ) dr^ dt^ 

<?x k k 

The components of the turbulent heat flow vector can be defined by 

c-p- u| T 1 = / u](F, Odt/uj^F htbwfF, F 1 , t. 1 1 ) d F^dt 1 


Since W is dependent on the difference r = t— t^, 

u | t1 = S' r d r lim “FT / u i (F .Odtfuj^F 1 , t- f)W(F, F t, t 1 ) dr 1 (r)dr = aX-^- 

dx k 0 2v . v 1 K dXLn »K lkj x . 


K Ik ^ = / u i ( r « l ) u k ( r 1 ,t - r) W (r , r hr) dr 1 

T OO rp 

The quantity a ( r ) = / K (r , r)dr means the eddy diffusivity of heat, 
o 

The concrete form of Green's functions V and W is determined by a geometry considered. 
In the simplest case when thedependence of the coefficients Eqs (8) and (20) upon the 
coordinates can be assumed negligible. Green's function may be found through the use of 
Fourier transformation. 

Then Green’s function for the velocity is given by 


I r — r 1 — U(t — t 1 ) 1 2 


v(f - r 1 , t-t 1 ) = 


4(i/+D)(t— t 1 ) 


2 \f 1 j[v + D)(t — t^P 


Green s function for the temperature W is the same, but we must use respectively a ^ 
instead of v, D. 
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The a approximation described above can be used to establish some semi-empirical methods 
of the calculation of the turbulent flow characteristics. However, in practice it is convenient 
to carry out the calculations using a more rough scheme. 

Consider fully developed flow in a circular tube. The equation for liquid motion is 


1 d , ,dU R 2 dp „rw n 
~r 77 f(^ + V JT = — • = - 2 — R 

f df 1 df p dz p 


(25) 


With the eddy diffusivity considered as a known function of the coordinates, after double 
integration we obtain 


r w R 1 A d A 

u(^)=— /— 4 - 

p v ? i 


(26) 


Using the common notations 


V2 _IL V * R - e 

* " p ’ ^ - ,0 ’ f ~ Hr 


we define 


^ 0 = 


U(€) 


% / 


1 fdf 


l + !l 


i 

/ . Vr \ 

JJlT 

0 l*3de 

o 1 + _T 


Re = 


2 U R n 2 } ^ 
— — = 2l £ / 


0 A i U T 

O l + I 


The friction factor £ is 


C—i 1 - F (Re) 

PL- ,JJ,2 

v * 


8 


V T . 

The function — = f(£ ( ^ Q ) being known, all flow characteristics may be calculated. A good 
expression for has been proposed by Reihardl [31 
326 
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— = 0.4 (r? - 11 th-I), 

V II 




According to this formula the eddy diffusivity of momentum Uj varies in the viscous layer 


as y d and for 77 >> 11 — ^0.477. Since in the turbulent core the velocity profile is 

v v 


described by the relation 


= 2.5 In 7 + 5.5 


the eddy diffusivity is to vary in the flow core as 


1 -± 


— ^ -0.4,f 

d<£ 


So it is reasonably to correct Reihardt’s relation and to determine the ratio of eddy and 
molecular viscosity as 


-j- = 0.4 [ v - 11 th^ ] ^ = 0.4 [t 7 0 (1 - f) - 11 th f 


The results of the velocity profile calculations based on Eqs (27) and (32) are given in Fig.l. 

In Fig. 2 the friction factor £ is plotted as a function of the number Re. 

The stabilized heat transfer coefficient calculation reduces to computation of the well-known 
Lyon integral 

1 , £ 9 


Nu =2! 0 


W ( /W(£)£d £) 2 


a prjp v ^ 1/ 

Assuming the turbulent Prandtl number equals to unity and using relation (32), the Nusselt 
number has been calculated in the range 0 < Pr 100, 5.10^ < Re < 10®. The results of the 
calculations are presented in Figs 3 and 4. 

In the figures these results are compared with data from the existing approximate equations. 
It is seen from this comparison that the approximation of the turbulent viscosity by relation (32) 
is satisfactory. 
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3. LOCAL HEAT TRANSFER IN A REACTOR CHANNEL 

As known, the heat transfer coefficient is dependent poorly upon a kind of the boundary 
conditions. However, in the nuclear reactors, where heat generation changes essentially 
it is necessary to take into account the effect of the axial heat flux variations on the heat 
transfer coefficient. 

Consider the heat transfer process in a circular tube. Let the heat flux q(z) change arbitrarily 
along the channel. What heat flux variation rate can result in local heat transfer coefficient 
and local temperature difference changes? 

Consider at first the ideal case of the 5-wise heat flux, i.e., the case when heat is gene- 
rated in an infinitely small element of thechannel length. How the temperature profile will 
be changed down the coolant flow? 

Obviously, in the heat source influence area the liquid near-wall layers temperature is 
much more than the flow core temperature. However, in time the liquid layers are mixing in 
the channel and the temperature profile is flattening. The length a: which temperature profile 
flattening occurs characterizes the temperature profile development rate. Let the function 
describing the temperature difference change along the channel is a heat flux influence func- 
tion and note it as G (z). 

Define the local temperature difference as 


(F w - T )cy Us * — f q (z^)G(z — zl)dz^ 

Z -oo 


T w — T x j . , 

r- r— = / f (x 1 ) G (x — x*) d x^ 


where 


ffx 1 ) = 2iiJ 
q(x) ’ 


If the heat flux is 5-wise, i.e., f(x) = S(x), then 


( T w _ T) 


cy U = G(x) 


OO 

The function G(x) can be normalized so that / G(x)dx = A. It is characterized by the 
1 » o _ 

relaxation length L = — / xG(x)dx and also by the high space moments x n =— f G(x) dx 
A o A o 

The influence function moments can be found by Fourier transformation. 
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Let 


„2 — 


G (p) = / e” px G (x)dx- / G (x)dx - p /xG(x)dx 4~— / x^G(x) dx - ... - A[l-pL + -- - ***] 

0 ° o (36) 


Then 


A = / G (x)dx *= G (p) I p-o 
o 


L = 


1 d G (p) . 

A dp P"*0’ 


1 i- j25 , 


(37) 


If a function f(x) changes poorly at distances of theorderof L, it can be expanded into a 
Taylor series and then 


T -T 
1 w 1 


St(x) 


, f f(x^)G(x - x^dx 1 =/f(x-x 1 )G(x 1 )dx 1 -f(x)7 G(x)dx-— / xG(x)dx + 

q(x) -oo o o dx o 


C y U 


(38) 


A[f(x)-L^ +...]- Af(x - L) 

dx 


Thus the local heat transfer coefficient is determined by the heat flux at the distance L up 
the flow. 

As we will see later, the factor A is the inverse Stanton number after full development 
of the temperature profile. 

In this connection 


t(x - L) 


(39) 


Thus using the stabilized heat transfer coefficient, we underestimate the heat transfer 
coefficient value when heat flux increases and overestimate it when heat flux decreases along 
a channel*. The last circumstance is essential for the water cooled reactors in which subcoolcd 
boiling is possible at core outlet. Therefore the fuel element surface temperature estimations 
and analysis of the subcooled boiling possibility performed in accordance with the usual heat 
transfer relations may be appealed quite optimistic. 

Surely, correctness of such estimations depends upon the rate of heat f|u^ change along a 
channel. If in expansion (38) the term Lrp can be neglected, i.e., L — ^ « 1, then the 

heat flux nonuniformity effect on heat transfer coefficient is small. Otherwise, if relative heat 
flux variations are small at distances about the relaxation length L, then this effect can be 
ignored. 

It is required to know the influence funcdon G (x) to calculate the local heat transfer coef- 
ficient for varying heat flux along a channel. The influence function G(x) can be determined 
from the solution of the energy equation. 


*It may be shown such a situation takes place in transient processes also, when q~const*f(t) 

326 
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This equation for liquid flowing in a circular tube is 


£ d £* a df a dx 2 dx 


and the boundary condition is 


!«r 

d£ 1 A 


( 10 ) 


(41) 


The solution Eq. (40) can be determined by means of the finite integral transformation over 
the eigenfunctions of SturnwLiouville’s equation 


with the boundary condition 

d Y. 


Id ,. a T d Y 


dY 


Tf f<1 + T' nr 


(42) 


df 


0. With the more general boundary condition 


^ ^ l| = *Y n (1) more complex problems may be solved. 

i 

Using the orthogonality of the function Y n , i.e., the equalities / £ w (<£) Y n Y m d<£ « 0 with 

o 

n ^ m, any function given in the interval [0.11 may be expanded into the series 


0 (f,x) « 1 


Y n (0 


1 


n /£w(f)Y 2 df 0 


/ £w(f)0(f,x)Y n (f) df 


(43) 


If to denote 


then the equality 


0 n (O * / fw (f) 0 (£, x) Y n (f)df, 
o 


(44) 


0(£,x) = S 


Y n (^)0 n (x) 


T 
n 1 


( 45 ) 


/£w(£)Y 2 d£ 


may be considered as the conversing formula for the finite integral transformation. 

Applying the integral transformation in the form of (44) to Eq. (40), i.e. multiplying the 
equation by the function Y n (f) and then integrating it from zero to unity, we obtain the equation 

_Pe _ ... q (x) R 

2 dx - n ' 


= Y n (l) — r— - /3 2 T n (x), 


(46) 


the solution of which is 


326 


T n (x) = Y 


O Q 1 

2R x . -jt/^U-x 1 ) i on 1 

„ U) _ ; q(xi)e Pe 


(47) 
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In accordance with Kq. (45) the temperature profile in the channel is determined by the series 

1 . Y n (f)Y n (l) x . -Itflx-x 1 ) 

I (f,x) = — 1 -j /q(x^)e Pe dx^ (4 

0 V° n n /fw(f)Y?df ‘°° 

o u 

Near the surface the temperature difference is 

T r 1 v Y n( X) X , -#^( X _ X 1) 

T W - ! “ ~ ff / q(x 1 ) e Pe dx 1 , (4! 

C>,lJ " /0 JfwWVjdf - 

o n 


and the inverse Stanton number 

1 T w -T x 
St(x) " q(x) 


Y 2 n (l) 

ts 




JtlxlMx 1 (50) 


Thus the heat flux influence function introduced before is 


Y 2 (1) 

G (x) » 2 JL_ e 

n ^°/ew(e)Y2d^ 

o u 

If the heat flux is constant, then the local Stanton number is 


__2 fl 2 
Pe^ n 


S Y n (1) 

St(x) 2 n/o o 1 0 

P n/fw(tfY 2 d£ 
o 


-|/9 2 x 
11 - e Pc " ] 


With x -* oo heat transfer is stabilized and 


1 Pe 


Y 2 (1) 


St(x) x- = St Q 2 n ^0 ~A I 

/9 2 /fw ( f)Y 2 df 

Calculating the moments of the function G(x) from Eq. (37) it can be shown that 


A = / G (x) dx — — 1 
o 2 d^O 


-11 - 


n <» 


Y 2 df 
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Yff(l) 


L “ St 0 7 xG(x)dx = ^ 
o ^ 


s '"'°/5 n 4 }fw({)v2a f 

' . y . 2 t‘> 


Limiting expansion (36) by the first two tenns, i.e., 


~ 1 11 
G (p )sf (1 - PL)— • — 

Stn St_ 1+pb 


we have the approximate representation of the function G(x) in the form 

G(x) - — 1 1. e (5 

St Q L 

It permits to calculate the local Stanton number for arbitraty heat flux distribution through 
the use of the relation 

~r pr 'cr r / f ( xl) e L (lxl (5 

St(x) St Q L _ M 


In particular for q - const in the thermal stabilization region 


(1- e 


The data on the heat transfer in the entrance section of tubes should be presented just in such a 

form. The thermal stabilization length may be considered as trebled relaxation length. It is 
S t 

corresponding to that — differs from unity by 5 per cent 
S t Q 

For practical calculations of the local heat transfer coefficient in case of complex heat 
flux distribution it is required to know the system of the eigenfunctions Y n (£). The calcula- 
tions such the functions and their eigenvalues have been performed usingthe code developed 
by A. I. Kuleshov and E.D.Beljaeva for computer M-20. The eddy and molecular diffusivity of 

heat ratio — cs Pr • — have been calculated using Eq. (32), and the velocity profile w(f) 
a v 

have been calculated using Eq. (27). The results of the thermal stabilization length calculations 
are given in Table 1. 


-12 - 
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Thermal stabilization length values 



P 

r = 0.01 


Pr = ( 

3.025 

Re 

i 

1 

10 5 

10 4 

10 5 


2.43 

8.5 

12.8 

5.01 

14.2 


JL, 

2 j 

Table 1 


Pr = 1.0 


In Fig. 5 Subbotin’s experimental data [7] on heat transfer for liquid metals in the entrance 
section of tubes are compared with the results of the calculations. The thermal stabilization 
length obtained from the experiments is in good agreement with the calculation. 

The local Nusselt number is presented in Fig. 6 for the sinusoidal heat flux. It is seen that 
for the decreasingpart of the sinusoid the number Nu is 20—40% less than the stabilized value, 

4. FRICTION AND HEAT TRANSFER IN PARALLEL FLOW THROUGH BUNDLES OF RODS 

Fuel assemblies of nuclear reactors are often made as bundles of cylindrical rods. Calcula- 
tions of friction and heat transfer in such channels are usually performed with the use of the 
relations obtained for flow in tubes. It is considered that the equivalent diameter defined as 
4S . 

djj = -pj- is adequate to account fora channel geometry. However, it has been shown in some 
experimental studies [8] that the equivalent diameter is not the universal parameter and does 
not provide simple describingof experimental data for bundles of rods. 

Therefore it is reasonable to carry out theoretical calculations of friction and heat transfer 
in such a geometry. 

For turbulent flow the full solution of the motion and energy equations can be obtained only 
numerically in such compoex geometry. Then some difficulties arise associated with the turbu- 
lent flow characteristics uncertainty. Therefore it is advisable to perform the calculation ap- 
proximately basing, however, on a reasonable description of turbulent flow properties. 

By such an approximate calculational method may be well-known in neutron physics Wigner- 
Seitz method in which a real lattice cell is replaced by some symmetric cylindrical one. It is 
obvious that such an operation can be applied only to extend bundles when it may be neglected 
by the azimuthal variations of shear stresses and the temperature difference over a rod perimeter. 
In such a symmetrized geometry heat transfer calculations have been made for liquid metal flow 
by Dwyer and Tee and also by Fridlander and Bonilla [9], 

a) Friction factor of a bundle of rods. 

Consider coolant flow in an annulus formed between the rod surface and the cell boundary. 

After integrating twice of the motion equation (25) with using of the boundary conditions we 
will have the velocity distribution over the annulus of the form. 

R 2 . ^ fdja-f 2 ) 

U,£ '--27 dT L ... -T (59) 


it f 

dz ii) 


wher £o ■ 
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The mean velocity over the annulus 


n 2 


U- - 


dp 

dz 1 (1 - £ 2 ) 2 d£ 


/ 


2/i(l— f 2 ) io £ ( 1 + — “ ) 


( 60 ) 


The eddy viscosity given — the velocity profile and friction factor can be calculated. 


r w 


D-2L 


-8 


(61) 


where 


4> - 


U 




$0 


1 (1 — £ 2 ) 2 d £ 


v* “ a-fo } (l-i 2 ) 2 $0 £(i+ — ) 

V 

The quantity rj Q = * ^ — is the characteristic parameter for turbulent flow in the 

annulus. The Reynolds number calculated through the equivalent diameter is determined by 
the relation 


Re 


U_d_ h 

v 


2>?o 


1 (l-^ 2 ) 2 ^ 
(l-£ 2 )(WJ 2 i f(l+ — ) 


(62) 


For laminar flow along a bundle of rods it can be made exact friction calculations. But these 
calculations are very complicated. Using Eqs (60), (61) and assuming i/'j' ■= 0 


dp 

dz 


U 


R z . e 3 

4/i 2 2 


lnr 

'IT' 


yU 2 

2g 


gjt 8 2 2 1-f 


yU 2 

2 Sdh 


(63) 


where 




64 


(_!_ ± ln * ) 
(1 — <)2 2 2 " 1 -e ’ 


Uy 

Mg 


X = 


dh x 


(64) 


2 

f m rnL m€ 2 

„R 2 0 

Let the quantity djj • x » d e be the effective diameter. With such a characteristic diameter 
definition the bundle friction law coincides with friction law for tubes, i.e. 
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64 _ 64 

He' x He e 


( 65 ) 


When f - 1, x -* -g- and £ -» as fora plane flat duct. 

It is required to know the eddy viscosity variation law over a channel section to calculate 
the friction factor of a bundle of rods. It may De assumed that the eddy viscosity distribution 
for a bundle of rods is described by the relation of the form 




- 11 th 


Vo (£-&)' 

1K1 -f 0 )_ 


1 -£ 

Wo 


( 66 ) 


Using relations (61) and (63) itmay be calculated the friction factor (as a function of the 
Reynolds number and the lattice density <. The results of the calculations are shown in Fig. 7. 

It may be seen that the coefficient £ is not a simple function of the number He, as it has 
been obtained from experiments. 

If the Reynolds number is defined through the effective diameter 


Re e = Re • x (67) 

then the exfoliation of the dependence £ = f (Re) with respect to the parameter * decreases 
essentially (1'ig. 7). These results permit to conclude that the effective diameter is a better 
geometrical characteristic of a bundle of rods than the equivalent diameter. However, it is true 
of the wide bundles only. 

b) The heat transfer coefficient for parallel flow through a bundle of rods. 

The equivalent cell method may be also used to calculate heat transfer in bundles of rods. 
The energy equation (40) integrating with the corresponding boundary conditions over the 
annulus, the expression for the Nusselt number will be found 


1 

Nu 


(l -$ 3 £ 0 £(l + Prl T ) 


( / fw(£)d £) 2 
£ 


( 68 ) 


Heat transfer calculations for a bundle of rods are performed using the velocity profile 
determined Eq. (59) and the eddy viscosity Eq. (66). 

The results of the calculations are presented in Figs 8,9,10. The Nusselt number calculated 
through the equivalent diameter is not a simple function of the Reynolds number as is the fric- 
tion factor £ 

If the effective diameter is introduced, then the exfoliation of the dependences 


Nu.f(Re), St = f (Re) 

decreases essentially. For this the data calculated are scattered by 10-15% from the criteria] 
relationship obtained for ducts. This result shows also that processing of data on heat transfer 
in bundles of rods should be made using the effective diameter and not the equivalent diameter. 
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NOMENCLATURE 
molecular diffusivity of heat 
pressure 

heat transfer perimeter 
tube radius, cell radius 
flow velocity 
temperature 

cross-sectional flow area 


u(r, t) = U + 

T(r, t) = T + t^ 
S 


x,y,z —coordinates 

y — specific weight 

S( 7) — delta-function 

p — viscosity 

p — liquid density, radius of rod 

v — kinematic viscosity 

f — dimensionless radial coordinate 

6=(-£)2 — ’’density" of lattice of rods 

R 

v* y 

n = — dimensionless distance from the wall 

v 

£ — friction factor 

r —shear stress 

Nu, Re, St, Pe, Pr — Nusselt, Reynolds, Stanton, Pecle, Prandtl numbers 




INDEXES 


w - wall 
T - turbulent 

0 - stabilized value 

1 - pulsational component 

- (dash) - notation of averaging 
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FIG. 3- FUNCTION Nu = f(Re, Pr) WITH Pr > 1 FOR 
LIQUID FLOW IN TUBES 

1 - Pt - 100; 2 - Pr = 10; 3 - Pr - 1.0; x Nu - 0.023 
Re 0.8 PjO.4 [4] ; y Nu - 0.021 Re 0 8 P r O-43 [5] 

Re ' Pr w\ 

o Nu = 16 J 

4.5(Pr^-l) , , • 

8 [ + 1.07] (1 82lg Re - 1.64) 2 

1.82 lg Re -1.64 




.-LL— 

-fH- 


Pi 

FIG.4. FUNCTION Nu - f(Pe) FOR LIQUID METAL FLOW IN 

TUBES 

1 - Nu » 7 + 0.025 Pe 0 ' 8 ; 2. Nu - 5-5 + 0.025 Pe 0 - 8 ; V Pr-0.005; 
x Pr - 0.01; o Pr - 0.025; A Pr - 0.05 

- 18 '• 
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FIG.7. FUNCTIONS £=f(Rc) AND ( = f(Re c )FOR A BUNDLE 


OF RODS 


1 . 0 -£ - 0.1 

2. -o- - f = 0.2 
$ — r = 0.3 


3. 't*, - f = 0.4 

^ - f = 0.5 

4. -p — ( = 0.6 


(1.82 lg Rc — 1.64) 2 
Points at left relate to Re e = Re*/ 

2i j € 3 Ja< 


p - rod radius, R — cell radius 
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FIG. 5. COMPARISON OF CALCULATIONAL AND EXPERIMENTAL DATA 
ON HEAT TRANSFER FOR LIQUID METAL FLOW IN THE ENTRANCE SEC- 
TION OF TUBES 

Experiment [7] o Pe = 500; A Pe = 1000; Calculation:! — Re = 9.67. 10 5 
Pr = 0.01; 2 - Re = 9-67.10 3 Pr = 0.025; 3 - Re =5 • I0 4 Pr = 0.01; 

4-Re = 10 5 Pr= 0.01; 5 -Re-10 5 ; Pr- 0.025; 



FIG.6. AS A FUNCTION OF RELATIVE CHANNEL LENGTH FOR 

Nu 0 

4 H 

SINUSOIDAL HEAT FLUX. Re = 10 , Pr = 1.0, - 50 

d 
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FIG. 8. HEAT TRANSFER CALCULATION FOR PARALLEL 
FLOW THROUGH A BUNDLE OF RODS FOR Pr = 1.0 


l.o - f = 0.1 

2.-0- - f = 0.2 

/f - ( = 0.3 
5. St = 


3. <*. - < = 0.4 

H - t = 0.5 

4. jp - ( = 0.6 

* - f = 0.7 


.7 V^(Pr^— 1)+1.07 


(1.82lgRe -1.64) 2 
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3. « = 0.4 


5. Nu = 7 + 0.025 Pe 0 ' 8 


Numbers Nu and Pe are calculated through equivalent diameter 

I Vt* of 
SOI 



FIG. 10. FUNCTION Nu e = f(Re e ) FOR PARALLEL LIQUID 
METAL FLOW OVER A BUNDLE OF RODS 

-o- - Pr = 0.005 o - P r = 0.01 A - Pr = 0.025 
x - Pr = 0.05 1. c = 0.1 2. t = 0.2 3. f = 0.4 

4. f = 0.6 5. Nu e = 7 + 0.025 Pe e °- 8 Nu e = Nu • f 
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Pe e = Pe • f 
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